/*
 * Zhihui Liu, last modify 20180202
 * init from Ethan Brodsky October 2011
**/

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <unistd.h>

#include "kmeans.h"
#define sqr(x) ((x)*(x))


#define MAX_ITERATIONS 100

#define KMEANS_DBG_ERR 1
#define KMEANS_DBG_ALERT 2
#define KMEANS_DBG_WARNING 3
#define KMEANS_DBG_INFO 4
#define KMEANS_DBG_TRACE 5
#define KMEANS_DBG_DEBUG 6
static int debug_level = KMEANS_DBG_INFO;
#define KMEANS_DBG(flag, level, fmt, args...)	\
{	\
	if (level <= debug_level )	\
	{	if (flag)	\
			printf("[kmeans dbg]<%s:%d> "fmt, __FILE__, __LINE__, ##args);	\
		else	\
			printf(fmt, ##args);	\
	}	\
}

#define KMEANS_ERR(fmt, args...)	\
	printf("[kmeans err]<%s:%d> "fmt, __FILE__, __LINE__, ##args)

#define STATIC_SYM

STATIC_SYM double km_calc_distance(int dim, double *p1, double *p2)
{
	double distance_sq_sum = 0;
	int i = 0;
	for (i = 0; i < dim; i++)
		distance_sq_sum += sqr(p1[i] - p2[i]);

	return distance_sq_sum;

}

STATIC_SYM double km_calc_abs_distance(int dim, double *p1, double *p2)
{
	double distance_sq_sum = 0;
	int i = 0;
	for (i = 0; i < dim; i++)
		distance_sq_sum += sqr(p1[i] - p2[i]);

	return sqrt(distance_sq_sum);

}


STATIC_SYM void km_calc_all_distances(int dim, int n, int k, double *X, double *centroid, double *distance_output)
{
	int i, j;
	for (i = 0; i < n; i++) {	
		for (j = 0; j < k; j++){ 
			//distance between point and cluster centroid
			distance_output[i*k + j] = km_calc_distance(dim, &X[i*dim], &centroid[j*dim]);
		}
	}
}

double km_rand(double n)
{
	float a = rand() / (RAND_MAX + 1.0);
	return n * a;
}
 
inline int
km_nearest(int dim, int n_cluster, double *x, double *centers, double *d2)
{
	int i;
	double d, min_d;
 
	for(i = 0; i < n_cluster; i++) {
		d = km_calc_distance(dim, x, &centers[i*dim]);
		if (i == 0) 
			min_d = d;
		if (min_d > d)
			min_d = d;
	}

	if (d2) *d2 = min_d;
	return 0;
}
 
int __km_centers_distances(int dim, int n, int n_cluster, double *X, double *centers, double *new_sum)
{
	double sum, *d = malloc(sizeof(double) * n);
	double min_d = 0;
	int j;
 
 	if (!d) {
		KMEANS_DBG(1, KMEANS_DBG_ERR, "malloc d failed\n");
		return -1;
	}

	sum = 0;
	/* compute the smallest dist for each node */
	for (j = 0; j < n; j++) {
		km_nearest(dim, n_cluster, &X[j*dim], centers, &d[j]);
		sum += d[j];
	}
	
	*new_sum = sum;
	/* choose the random node */
	sum = km_rand(sum);
	for (j = 0; j < n; j++) {
		if ((sum -= d[j]) > 0) 
			continue;
		break;
	}

	free(d);
	return j;
}

int km_init_centers(int dim, int n, int n_cluster, double *X, double *centers)
{
	int i, j, m, ret = 0;
	int rand_index = 0, index = 0;
	double sum, min_sum;
	double min_d = 0;
	double min_centers[MAX_DIM];
	int trial = 2 + (int)log(n_cluster);
    //trial = 10;
 
	srand(time(NULL));
 	rand_index = (int)km_rand(n);
	KMEANS_DBG(1, KMEANS_DBG_INFO, "init sample(cnt:%d, dim:%d) center 0 index is %d\n", n, dim, rand_index);
 	for (i = 0; i < dim; i++) {
		centers[i] = X[rand_index *dim + i];
	}

	for (i = 1; i < n_cluster; i++) {
		sum = 0;
		/* compute the smallest dist for each node */
		for (j = 0; j < trial; j++) {
			index = __km_centers_distances(dim, n, i, X, centers, &sum); 
			if (index < 0) {
				KMEANS_DBG(1, KMEANS_DBG_ERR, "init center failed\n");
				ret = -1;
				goto out;
			}
			/* set the next center */
	 		for (m = 0; m < dim; m++) { 
				centers[i*dim + m] = X[index *dim + m];
				if (j == 0) {
					min_centers[m] = X[index *dim + m];
				}
			}

			/* compute dist after adding a new center */
			index = __km_centers_distances(dim, n, i+1, X, centers, &sum); 
			if (j == 0)
				min_sum = sum;
			if (sum < min_sum) {
				min_sum = sum;
		 		for (m = 0; m < dim; m++) 
					min_centers[m] = centers[i*dim + m];
			}

		}
		/* use min centers */
		for (m = 0; m < dim; m++) 
			centers[i*dim + m] = min_centers[m];

	}

	KMEANS_DBG(0, KMEANS_DBG_INFO, "[init centers] n_cluster:%d\n", n_cluster);
	for(i = 0; i < n_cluster; i++) {
		KMEANS_DBG(0, KMEANS_DBG_INFO, "center[%d]:",  i);
		for (j = 0; j < dim; j++) {
			KMEANS_DBG(0, KMEANS_DBG_INFO, " %16f", centers[i*dim +j]);
		}
		KMEANS_DBG(0, KMEANS_DBG_INFO, "\n");
	}

out:
	return ret;
}
 
// NOTE: a point with cluster assignment -1 is ignored
STATIC_SYM double km_calc_total_distance(int dim, int n, int k, double *X, double *centroids, int *cluster_assignment_index)
{
	double tot_D = 0;
	int i = 0, c_index = 0;

	for (i = 0; i < n; i++)
	{
		c_index = cluster_assignment_index[i];
		if (c_index == -1)
			continue;
		tot_D += km_calc_distance(dim, &X[i*dim], &centroids[c_index*dim]);
	}

	return tot_D;
}

void km_choose_all_clusters_from_distances(int dim, int n, int k, double *distance_array, int *X_cluster_index)
{
	int i, j;
	int index = -1;
	double min_dist, cur_dist;

	for (i = 0; i < n; i++) {
		for (j = 0; j < k; j++) {
			// distance between point and cluster centroid
			cur_dist = distance_array[i*k + j];
			if (j == 0) {
				min_dist = cur_dist;
				index = 0;
			}
			if (cur_dist < min_dist)
			{
				index = j;
				min_dist = cur_dist;
			}
		}

		// record in array
		X_cluster_index[i] = index;
	}
}

void km_calc_centers_k_mediod(int dim, int n, int k, double *X, int *X_cluster_index, double *new_centers)
{
	int cluster_member_count[MAX_CLUSTERS];
	int i, j, c_index;

	memset(cluster_member_count, 0x0, sizeof(int)*MAX_CLUSTERS);
	memset(new_centers, 0x0, sizeof(double)*k*dim);

	for (i = 0; i < n; i++) {
		c_index = X_cluster_index[i];
		cluster_member_count[c_index]++;

		//sum point coordinates for finding centroid
		for (j = 0; j < dim; j++)
			new_centers[c_index*dim + j] += X[i*dim + j];
	}


	//divide coordinate sum by number of members to find centroid
	for (i = 0; i < k; i++) {
		if (cluster_member_count[i] == 0) {
			KMEANS_DBG(1, KMEANS_DBG_ERR, "WARNING: Empty cluster %d! \n", i);
			continue;
		}

		for (j = 0; j < dim; j++)
			new_centers[i*dim + j] /= cluster_member_count[i];  /// XXXX will divide by zero here for any empty clusters!

	}
}


void km_calc_centers(int dim, int n, int k, double *X, int *X_cluster_index, double *new_centers)
{
	int cluster_member_count[MAX_CLUSTERS];
	int i, j, c_index;

	memset(cluster_member_count, 0x0, sizeof(int)*MAX_CLUSTERS);
	memset(new_centers, 0x0, sizeof(double)*k*dim);

	for (i = 0; i < n; i++) {
		c_index = X_cluster_index[i];
		cluster_member_count[c_index]++;

		//sum point coordinates for finding centroid
		for (j = 0; j < dim; j++)
			new_centers[c_index*dim + j] += X[i*dim + j];
	}


	//divide coordinate sum by number of members to find centroid
	for (i = 0; i < k; i++) {
		if (cluster_member_count[i] == 0) {
			KMEANS_DBG(1, KMEANS_DBG_ERR, "WARNING: Empty cluster %d! \n", i);
			continue;
		}

		for (j = 0; j < dim; j++)
			new_centers[i*dim + j] /= cluster_member_count[i];  /// XXXX will divide by zero here for any empty clusters!

	}
}

void km_get_cluster_member_count(int n, int k, int *labels, int *cluster_member_count)
{
	int i;
	memset(cluster_member_count, 0x0, sizeof(int) *k);
	for (i = 0; i < n; i++)
		cluster_member_count[labels[i]]++;
}

void km_cluster_diag(int dim, int n, int k, double *X, int *cluster_assignment_index, double *centers)
{
	int cluster_member_count[MAX_CLUSTERS];
	int i, j;
	char buf[256];
	int offset = 0;

	km_get_cluster_member_count(n, k, cluster_assignment_index, cluster_member_count);
	KMEANS_DBG(1,KMEANS_DBG_INFO,"Final clusters:\n");
	for (i = 0; i < k; i++) {
		offset = 0;
		offset += snprintf(buf, sizeof(buf), "\tcluster %d:     members count: %8d, centroid (", i, cluster_member_count[i]);
		for (j = 0; j < dim; j++) {
			offset += snprintf(buf + offset, sizeof(buf)-offset, "%16f ", centers[i*dim + j]);
		}
		snprintf(buf + offset, sizeof(buf)-offset, ")");
		KMEANS_DBG(1, KMEANS_DBG_INFO,"%s\n", buf);
	}		
}

void km_copy_assignment_array(int n, int *src, int *dst)
{
	memcpy(dst, src, sizeof(int)*n);
}

int km_assignment_change_count(int n, int a[], int b[])
{
	int change_count = 0;
	int i = 0;

	for (i = 0; i < n; i++)
		if (a[i] != b[i])
			change_count++;

	return change_count;
}


void km_quicksort(double a[], int left, int right)
{
	int i, j;
	double x;
	if (left < right) {
		i = left, j = right, x = a[left];
		while (i < j) {
			while(i < j && a[j] >= x)
				j--;  
			if(i < j) 
				a[i++] = a[j];

			while(i < j && a[i] < x) 
				i++;  
			if(i < j) 
				a[j--] = a[i];
		}
		a[i] = x;
		km_quicksort(a, left, i - 1); 
		km_quicksort(a, i + 1, right);
	}
}

int _km_percentile(double *X, int n, double percent, double *q )
{
        double index = 0;
        unsigned int lo = 0;
        double w1 = 0,w2 = 0,hi = 0;

        index = ( n - 1 ) * (double)percent/100 + 1;
        if (index <= 1) {
                *q = X[0];
                goto out;
        }

        if (index >= n) {
                *q = X[n-1];
                goto out;
        }

        lo = (unsigned int)index;
        hi = index - (double)lo;
        lo = lo - 1;

        w1 = X[lo];
        w2 = X[lo+1];

        *q = w2 * hi + w1 * ( 1 - hi );
out:
	return 0;
}


int km_test_data(int n, double *X)
{
	FILE *fp = NULL;
	fp = fopen("percentile_test.txt", "w");
	int i = 0;
	for (i =0; i < n; i++)
		fprintf(fp, "%f 123456\n", X[i]);

	fclose(fp);
	return 0;
}

int km_percentile(int dim, int n, int k, double *X, kmeans_model_t *m, int *labels)
{
	int i, j;
	//double *centers = m->centers;
    kmeans_cluster_t *cluster = NULL;
	struct kmeans_percentile *d = NULL;
	double *dists;
	int cnt;

	dists = malloc(sizeof(double)*n);
	if (!dists) {
		printf("malloc dists failed\n");
		return -1;
	}

	for (i = 0; i < k; i++) {
		cnt = 0;
        cluster = &m->clusters[i];
        d = &cluster->dist;
		for(j = 0; j < n; j++) {
			if (labels[j] == i) {
				dists[cnt] = km_calc_abs_distance(dim, &X[j*dim], cluster->center);
				cnt++;
			}
		}
		if (cnt <= 0) {
			KMEANS_DBG(1, KMEANS_DBG_ALERT, "[ALERT]:cluster %d has no sample\n", i);
			d->q1 = d->q2 = d->q3 = d->q4 = d->q5 = d->q6 = d->q7 =0;
            for (j = 0; j <= 10; j++) {
                d->quantile[j] = 0;
            }
			KMEANS_DBG(1, KMEANS_DBG_INFO, "cluster[%d],cnt:%d, q1:%f,q2:%f,q3:%f,q4:%f,q5:%f,q6:%f,q7:%f\n", 
				i,cnt,  d->q1, d->q2, d->q3, d->q4, d->q5, d->q5, d->q7);

			continue;
		}
		km_quicksort(dists, 0, cnt - 1);
		//km_test_data(cnt, dists);
		_km_percentile(dists, cnt, 25, &d->q1);
		_km_percentile(dists, cnt, 50, &d->q2);
		_km_percentile(dists, cnt, 75, &d->q3);
		d->q4 = d->q1 - 1.5 *(d->q3 - d->q1);
		d->q5 = d->q3 + 1.5 *(d->q3 - d->q1);
		d->q6 = d->q1 - 3 *(d->q3 - d->q1);
		d->q7 = d->q3 + 3 *(d->q3 - d->q1);

        for (j = 0; j <= 10; j++) {
            _km_percentile(dists, cnt, 100 - 0.1 * (10-j), &d->quantile[j]);
        }

		KMEANS_DBG(1, KMEANS_DBG_INFO, "cluster[%d],cnt:%d, q1:%f,q2:%f,q3:%f,q4:%f,q5:%f,q6:%f,q7:%f, quantile:(", 
			i,cnt,  d->q1, d->q2, d->q3, d->q4,d->q5,d->q5,d->q7);
        for (j = 0; j <= 10; j++) {
            KMEANS_DBG(0, KMEANS_DBG_INFO, "%f, ", d->quantile[j]);
        }
        KMEANS_DBG(0, KMEANS_DBG_INFO, ")\n");

	}

	if (dists)
		free(dists);
	return 0;
}

/*
 *	Return:
 *		1, has zero cluster;
 *		0, OK
 *		-1, error
 */
int kmeans(int dim, double *X, int n, int k, 
		kmeans_model_t *m, int *labels)
{
	int ret = 0, change_count = 0;
	double totD = 0, prev_totD = 0;
	int batch_iteration = 0;
	double centers[MAX_DIM * MAX_CLUSTERS];
	int cluster_count[MAX_CLUSTERS];
	int i, j;

	double *dist = (double *)malloc(sizeof(double) * n * k);
	int *labels_cur  = (int *)malloc(sizeof(int) * n);
	int *labels_prev = (int *)malloc(sizeof(int) * n);

	KMEANS_DBG(1, KMEANS_DBG_TRACE, "start kmeans for data dim %d to cluster %d\n", n, k);
	if (!dist || !labels_cur || !labels_prev) {
		KMEANS_ERR("Error allocating dist arrays, dist:%p, cur:%p, prev:%p\n", 
				dist, labels_cur, labels_prev);
		ret = -1;
		goto out;
	}

	m->dim = dim;
	m->n_cluster = k;

    memset(centers, 0x0, sizeof(centers));
    memset(cluster_count, 0x0, sizeof(cluster_count));
	// initial setup  
	if (km_init_centers(dim, n, k, X, centers)) {
		KMEANS_DBG(1, KMEANS_DBG_ERR, "init center failed\n");
		ret  = -1;
		goto out;
	}
	km_calc_all_distances(dim, n, k, X, centers, dist);
	km_choose_all_clusters_from_distances(dim, n, k, dist, labels_cur);
	km_copy_assignment_array(n, labels_cur, labels_prev);

	// BATCH UPDATE
	while (batch_iteration < MAX_ITERATIONS) {
		KMEANS_DBG(1, KMEANS_DBG_TRACE, "batch iteration %d \n", batch_iteration);
		km_calc_centers(dim, n, k, X, labels_cur, centers);

		//deal with empty clusters
		//XXXXXXXXXXXXXX

		//check converged
		totD = km_calc_total_distance(dim, n, k, X, centers, labels_cur);
		if (prev_totD == 0)
			prev_totD = totD;
		if (totD > prev_totD) {
			KMEANS_DBG(1, KMEANS_DBG_INFO, "negative progress made on this step - iteration completed (%.2f) \n", totD - prev_totD);
			km_copy_assignment_array(n, labels_prev, labels_cur);
			km_calc_centers(dim, n, k, X, labels_cur, centers);
			break;
		}

		// save previous step
		km_copy_assignment_array(n, labels_cur, labels_prev);

		// move all points to nearest cluster
		km_calc_all_distances(dim, n, k, X, centers, dist);
		km_choose_all_clusters_from_distances(dim, n, k, dist, labels_cur);

		change_count = km_assignment_change_count(n, labels_cur, labels_prev);
		KMEANS_DBG(1, KMEANS_DBG_TRACE, "inter times:%3d, cluster changed number:%9d, new cost:%16.2f, new cost - old cost:%17.2f\n", batch_iteration, change_count, totD, totD - prev_totD);

		// done with this phase if nothing has changed
		if (change_count == 0) {
			KMEANS_DBG(1, KMEANS_DBG_TRACE, "no change made on this step - iteration completed \n");
			break;
		}

		prev_totD = totD;
		batch_iteration++;
	}


	km_cluster_diag(dim, n, k, X, labels_cur, centers);
	//don't support ONLINE UPDATE
	/* The online update prtion of this code has never worked properly, but batch update has been adequate for our projects so far. */

	ret = 0;
	// write to output array
	km_copy_assignment_array(n, labels_cur, labels);    
	km_get_cluster_member_count(n, k, labels, cluster_count);
    for (i = 0; i < k; i++) {
        for (j = 0; j < dim; j++) 
            m->clusters[i].center[j] = centers[i*dim + j];
        m->clusters[i].sample_count = cluster_count[i];
    }

	for (i = 0; i < k; i++) {
		if (cluster_count[i] == 0) {
			KMEANS_DBG(1, KMEANS_DBG_WARNING, "cluster(total :%d) index %d with zero node\n", k, i);
			ret = 1;
		}
	}
	if (!ret)
		km_percentile(dim, n, k, X, m, labels);
out:
	if (dist)
		free(dist);
	if (labels_cur)
		free(labels_cur);
	if (labels_prev)
		free(labels_prev);

	return ret;
}           

int kmeans_check_preprocess(kmeans_model_t  *m, double *s)
{
    int j = 0;
    for(j = 0; j < m->dim; j++) {
        s[j] = s[j] * m->dims[j].enlarge_factor;
    }
    return 0;
}

/*
 * type = 0, quantile; type = 1, definate;
*/
int kmeans_check(kmeans_model_t  *m, int n, double *samples, double potential, double definate, int type)
{
	int anomaly_cnt = 0, i ,j= 0, k = 0, index;
	double min_dist = -1, dist = 0;
	double q4 = 0;
	int dim = m->dim;
	int ano_cnt = 0, definate_cnt = 0;
    double check_potential_dist = 0;
    double check_definate_dist = 0;
    kmeans_cluster_t *c = NULL;
    double s[MAX_DIM];

    if (!m || n <= 0 || !samples) {
        printf("Invalid args for kmeans check\n");
        return -1;
    }

	for (i = 0; i < n; i++) {
		min_dist = -1;
        memcpy(s, &samples[i*dim], sizeof(double) * dim);
        kmeans_check_preprocess(m, s);

        /* find cluster */
		for(j = 0; j < m->n_cluster; j++) {
            c = &m->clusters[j];
			if (c->sample_count <= 0) {
				/* ignore zero cluster */
				continue;
			}
			dist = km_calc_abs_distance(dim, c->center, s);
			if (min_dist < 0 || min_dist > dist) {
				min_dist = dist;
				index = j;
			}
		}
        
        c = &m->clusters[index];
        if (type == 0) {
            check_potential_dist = c->dist.quantile[(int)( (1- potential) * 10)];
            check_definate_dist = c->dist.quantile[(int)(1* 10 - definate * 10)];
        } else {
            check_potential_dist = c->dist.q5;
            check_definate_dist = c->dist.q7;
        }

		if  (min_dist > check_definate_dist) {
            printf("=============>def:%16f, index :%16f, minist %16f, check_definate_dist:%16f\n", 
                definate, ((1 - definate) * 10), min_dist, check_definate_dist);
            definate_cnt++;
			ano_cnt++;
			KMEANS_DBG(0, KMEANS_DBG_DEBUG, "Anomaly definate sample(n_cluster:%d,sample:%d,dist:%16f, checkdist:(%16f,%16f)), sample:", 
				index, i, min_dist, check_potential_dist, check_definate_dist);
			for (k = 0; k < dim; k++)
				KMEANS_DBG(0, KMEANS_DBG_DEBUG,"%16f", samples[i*dim + k]);
			KMEANS_DBG(0, KMEANS_DBG_DEBUG,"\n");
		} else if (min_dist > check_potential_dist) {
            ano_cnt++;
			KMEANS_DBG(0, KMEANS_DBG_DEBUG, "Anomaly potential sample(n_cluster:%d,sample:%d,dist:%16f, checkdist:(%16f,%16f)), sample:", 
				index, i, min_dist, check_potential_dist, check_definate_dist);
			for (k = 0; k < dim; k++)
				KMEANS_DBG(0, KMEANS_DBG_DEBUG,"%16f", samples[i*dim + k]);
			KMEANS_DBG(0, KMEANS_DBG_DEBUG,"\n");

        } else {
			KMEANS_DBG(0, KMEANS_DBG_DEBUG, "Nomaly sample(n_cluster:%d,sample:%d,dist:%16f, q5:%16f), sample:", 
				index, i, min_dist, c->dist.q5);
			for (k = 0; k < dim; k++)
				KMEANS_DBG(0, KMEANS_DBG_DEBUG,"%16f", samples[i*dim + k]);
			KMEANS_DBG(0, KMEANS_DBG_DEBUG,"\n");
		}
	}
	KMEANS_DBG(0, KMEANS_DBG_INFO,"[ CHECK ] total cnt:%d, ano_cnt:%d, def_cnt :%d, ano percent:%f, def_percert:%f\n", 
        n, ano_cnt, definate_cnt, (double)ano_cnt/n, (double)definate_cnt/n);
	return ano_cnt;
}

int __kmeans_silhouette(int dim, double *X, int n, kmeans_model_t *m, int *labels, double *cluster_si)
{
	double d;
	double *sis = NULL;
	double total_si = 0;
	double ai;
	double bi[MAX_CLUSTERS], min_bi;
	int bi_cnt[MAX_CLUSTERS];
	int i = 0, j = 0;
	int ret = 0;
	int first = 0;

#if 0
	dobule *dist;
	dist = malloc(sizeof(double) * n * n);
	if (!dist) {
		KMEANS_DBG(1, KMEANS_DBG_ERR, "malloc dist failed for silhouette\n");
		ret = -1;
		goto out;
	}
#endif
	sis = malloc(sizeof(double) *n);
	if (!sis) {
		KMEANS_DBG(1, KMEANS_DBG_ERR, "malloc silhouette failed for silhouette\n");
		ret = -1;
		goto out;
	}

	for (i = 0; i < n; i++ ) {
		/* compute all bi for i node */
		memset(bi, 0x0, sizeof(double) * MAX_CLUSTERS);
		memset(bi_cnt, 0x0, sizeof(int) * MAX_CLUSTERS);
		for (j = 0; j < n; j++) {
			d = km_calc_abs_distance(dim, &X[i*dim], &X[j*dim]);
			//dist[i *n + j] = d;
			bi[labels[j]] += d;
			bi_cnt[labels[j]]++;
		}

		/* get ai and smallest bi */
		for (j = 0; j < MAX_CLUSTERS; j++) {
			if (j == labels[i]) {
				bi_cnt[j]--;	/* minus i itself */
				ai = (bi_cnt[j] == 0) ? 0: bi[j]/bi_cnt[j];
				continue;
			}
			bi[j] = (bi_cnt[j] == 0) ? 0: bi[j]/bi_cnt[j];
		}

		min_bi = 0; first = 0;
		for (j = 0; j < m->n_cluster; j++ ) {
			if  (j == labels[i])
				continue;
			if (first == 0) {
				min_bi = bi[j];
				first = 1;
			}
			if (min_bi > bi[j])
				min_bi = bi[j];
		}

		/* compute si */
		if (ai == 0 && min_bi == 0) {
			KMEANS_DBG(1, KMEANS_DBG_WARNING, "cluster cnt %d node :%d silhouette is zero, ai:%f, min_bi:%f\n", m->n_cluster, i, ai, min_bi);
			sis[i] = 0;
		} else  {
			sis[i] =  (min_bi - ai) / (min_bi > ai? min_bi: ai);
		}
		total_si += sis[i];
	}

	*cluster_si = total_si / n;
	KMEANS_DBG(1, KMEANS_DBG_INFO, "[silhouette] si for cluster %d is :%f, total:%f\n", m->n_cluster, *cluster_si, total_si); 

out:
#if 0
	if (dist) {
		free(dist);
	}
#endif
	if (sis)
		free(sis);

	return ret;
}

int kmeans_silhouette(int dim, double *X, int n, int *best_cluster)
{
	int i, j;
	double silhouette[MAX_CLUSTERS];
	int cluster_count[MAX_CLUSTERS];
	int n_start = 2;
	kmeans_model_t m;
	int *labels;
	int ret = 0;
	double max_si = 0;
	int max_cluster = 0;
	int n_cluster = 0;

	labels = malloc(sizeof(int) * n);
	if (!labels) {
		KMEANS_DBG(1, KMEANS_DBG_ERR,"create labels failed for silhouette\n");
		ret = -1;
		goto out;
	}

	for (i = n_start; i <= MAX_CLUSTERS; i++) {
		memset(&m, 0x0, sizeof(kmeans_model_t));
		KMEANS_DBG(1, KMEANS_DBG_INFO,"===> kmeans compute silhouette for n_cluster %d\n", i);
		ret = kmeans(dim, X, n, i, &m, labels);
		if (ret == 1) {
			KMEANS_DBG(1, KMEANS_DBG_INFO, "max cluster is %d\n", max_cluster);
			break;
		}
		if (ret < 0) {
			KMEANS_DBG(1, KMEANS_DBG_ERR,"kmeans failed for n_cluster %d\n", i);
			ret = -1;
			goto out;
		}
		ret = __kmeans_silhouette(dim, X, n, &m, labels, &silhouette[i-1]);
		if (ret) {
			KMEANS_DBG(1, KMEANS_DBG_ERR,"kmeans compute silhouette failed for n_cluster %d\n", i);
			ret = -1;
			goto out;
		}
		KMEANS_DBG(1, KMEANS_DBG_INFO, "si for cluster %d is :%f\n", i, silhouette[i-1]); 
		max_cluster = i;
	}
	if (max_cluster < n_start) {
		KMEANS_DBG(1, KMEANS_DBG_ALERT,"kmeans compute only one cluster\n");
		*best_cluster = 1;;
		goto out;
	}
	for (i = n_start; i <= max_cluster; i++) {
		KMEANS_DBG(1, KMEANS_DBG_INFO,"n_cluster %d silhouette %f\n", i, silhouette[i-1]);
		if (i == n_start) {
			max_si = silhouette[i-1];
			n_cluster = i;
		}
		if (silhouette[i-1] > max_si) {
			max_si = silhouette[i-1];
			n_cluster = i;
		}
	}
	*best_cluster = n_cluster;
	KMEANS_DBG(1, KMEANS_DBG_ERR,"[best] n_cluster %d silhouette %f\n", *best_cluster, silhouette[*best_cluster-1]);

	ret = 0;
out:
	if (labels) {
		free(labels);
	}
	
	return ret;
}

int kmeans_sample_preprocess(int dim, double *X, int n, kmeans_model_t *m)
{ 
    int i = 0, j = 0;
    double means[MAX_DIM];
    double max_mean = 0, factor = 0;

    memset(means, 0x0, sizeof(double) * MAX_DIM);
    for (i = 0; i < n; i++) {
        for(j = 0; j < dim; j++) {
            means[j] += X[i*dim + j]; 
        }
    }

    for(j = 0; j < dim; j++) {
        means[j] = means[j]/n;
        m->dims[j].mean = means[j];
        
        if (max_mean == 0)
            max_mean = means[j];

        if (means[j] > max_mean)
            max_mean = means[j];
    }

    for(j = 0; j < dim; j++) {
        m->dims[j].enlarge_factor = 1;
        if (m->dims[j].mean == 0) {
            KMEANS_DBG(1, KMEANS_DBG_TRACE, "kmeans dim %d mean %f enlarge_factor %f\n", 
                j, m->dims[j].mean, m->dims[j].enlarge_factor);
            continue;
        }
        
        factor = max_mean/means[j];
        if (factor >= 5)
            m->dims[j].enlarge_factor = factor; 

        KMEANS_DBG(1, KMEANS_DBG_TRACE, "kmeans dim %d mean %f enlarge_factor %f\n", 
            j, m->dims[j].mean, m->dims[j].enlarge_factor);
    }

    for (i = 0; i < n; i++) {
        for(j = 0; j < dim; j++) {
            X[i*dim+j] = X[i*dim+j] * m->dims[j].enlarge_factor;
        }
    }

    return 0;
}

int kmeans_noisy()
{
    /* use DBSCAN to remove noisy node first */
    return 0;
}

int kmeans_train(int dim, double *X, int n, kmeans_model_t *m, int best_n_cluster)
{
	int n_cluster = 0;
	int *labels = NULL;
	int i, j;
	int *cluster_count;
	int ret = 0;
    kmeans_cluster_t *cluster = NULL;

	if (!X || n <= 0 || !m || dim <= 0) {
		KMEANS_DBG(1, KMEANS_DBG_ERR, "invalid arg\n");
		ret = -1;
		goto out;
	}

	KMEANS_DBG(1, KMEANS_DBG_INFO, "kmeans train start with best_n_cluster set to %d===>\n", best_n_cluster);	
	memset(m, 0x0, sizeof(kmeans_model_t));
    kmeans_sample_preprocess(dim, X, n, m);

	if (best_n_cluster)
		n_cluster = best_n_cluster;
	else  {
		ret = kmeans_silhouette(dim, X, n, &n_cluster);
		if (ret) {
			printf("get best cluster failed\n");
			ret = -1;
			goto out;
		}
		printf("best cluster is :%d\n", n_cluster);
	}

	if (n_cluster == 1) {
		KMEANS_DBG(1, KMEANS_DBG_ERR, "only one cluster create\n");	
		goto out;
	}

	labels = malloc(sizeof(int) * n);
	if (!labels) {
		printf("malloc labels failed\n");
		ret = -1;
		goto out;
	}
	
	ret = kmeans(dim, X, n, n_cluster, m, labels);
	if (ret == 1) {
		KMEANS_DBG(1, KMEANS_DBG_WARNING, "cluster with 0 node created\n");
	}
	if (ret < 0) {
		printf("kmeans failed\n");
		ret = -1;
		goto out;
	}

	j = 0;
	for (i = 0; i < n; i++) {
		KMEANS_DBG(0, KMEANS_DBG_DEBUG, "%d ", labels[i]);
		j++;
		if (j == 20) {
			KMEANS_DBG(0, KMEANS_DBG_DEBUG, "\n");
			j = 0;
		}
	}

	KMEANS_DBG(0, KMEANS_DBG_INFO, "[train sample model] dim:%d, n_cluster:%d\n", m->dim, m->n_cluster);
	for (j = 0; j < dim; j++) {
		KMEANS_DBG(0, KMEANS_DBG_INFO, "\tdim %d mean %16f enlarge_factor %16f\n", m->dims[j].mean, m->dims[j].enlarge_factor);
	}
	for(i = 0; i < n_cluster; i++) {
        cluster = &m->clusters[i];
		KMEANS_DBG(0, KMEANS_DBG_INFO, "cluser(cnt:%6d) %d center:", cluster->sample_count, i);
		for (j = 0; j < dim; j++) {
			KMEANS_DBG(0, KMEANS_DBG_INFO, " %16f", cluster->center[j]);
		}
		KMEANS_DBG(0, KMEANS_DBG_INFO, "\n");
		KMEANS_DBG(0, KMEANS_DBG_INFO,"\tq1:%16f, q2:%16f, q3:%16f, q4:%16f, q5:%16f, q6:%16f, q7:%16f, quantile(",
		    cluster->dist.q1, cluster->dist.q2, cluster->dist.q3,
			cluster->dist.q4, cluster->dist.q5,
			cluster->dist.q6, cluster->dist.q7);
        for (j = 0; j < sizeof(cluster->dist.quantile)/sizeof(cluster->dist.quantile[0]); j++) {
            KMEANS_DBG(0, KMEANS_DBG_INFO, "%f, ", cluster->dist.quantile[j]);
        }
        KMEANS_DBG(0, KMEANS_DBG_INFO, ")\n");


	}

out:

	if (labels) 
		free(labels);

	return ret;
}


#ifdef KMNEAS_TEST
#define MAX_TRAIN_SAMPLE 10000
double km_samples[MAX_TRAIN_SAMPLE] = {0};
int km_sample_cnt = 0;
int dim_2_is_string = 1;

int km_treate_sample(char *buf, int dim, double *sample)
{
	char *token;
	const char s[2] = " ";
	int len = strlen(buf);
	while(buf[len -1] == '\r' || buf[len -1] == '\n') {
		len--;
		if (len < 0) {
            printf("sample len 0\n");
			return -1;
		}
	}
	
	token = strstr(buf, s);
	if (!token) {
        token = strstr(buf, "\t");
        if (!token) {
            printf("no space or tab in sample\n");
		    return -1;
        }
    }
	*token = '\0';
	sample[0] = strtod(buf, NULL);
    //sample[0] = sample[0] * 10;

	if (dim == 1)
		return 0;
	token++;
	while ((*token == ' ' || *token == '\t')&& token < buf+len) {
		token++;
	}
	if (token == buf + len)
		return -1;

    if (dim_2_is_string)
	    sample[1] = strlen(token);
    else
        sample[1] = strtod(token, NULL); 

	return 0;
}

int km_read_samples(char *filename, int dim)
{
	FILE *fp;
	char buf[1024] = {0};
	int cnt = 0;
	char *p;
	double s[MAX_DIM];
	int i, j;
    double count[MAX_DIM];
	
	if (dim > 2) {
		printf("km_treate_sample only support <2 dim now\n");
		return -1;
	}

	if ((fp = fopen(filename, "r")) == NULL) {
		return -1;
	}
	
	while(1) {
		p = fgets(buf, sizeof(buf), fp);
		if(!p){
			break;
		}
		if (km_treate_sample(buf, dim, s)) {
			printf("treat sample failed, buf:%s\n", buf);
			continue;
		}
		for (i = 0; i < dim; i++) {
			km_samples[cnt * dim + i] = s[i];
		}
		cnt++;
		if (cnt >= MAX_TRAIN_SAMPLE)
			break;
	}

	km_sample_cnt = cnt;
    memset(count, 0x0, sizeof(count));
	for (i = 0; i < cnt; i++){
		KMEANS_DBG(0, KMEANS_DBG_DEBUG,"index: %d data:", i);
		for (j = 0; j < dim; j++) {
			KMEANS_DBG(0, KMEANS_DBG_DEBUG, " %16f", km_samples[i *dim + j]);
            count[j] += km_samples[i *dim + j];
		}
		KMEANS_DBG(0, KMEANS_DBG_DEBUG, "\n");
	}
	fclose(fp);
    for (j = 0; j< dim; j++) {
        KMEANS_DBG(1, KMEANS_DBG_INFO,"sample dim %d mean %f\n", j, count[j]/km_sample_cnt);
    }

	return cnt;
}


char check_filename[256] = {0,};
char train_filename[256] = {0,};
int test_dim = 1;
int best_n_cluster = 0;
int boxplot_check = 0;
int _parse_opt(int argc, char **argv)
{
        char c;
        //const char optstr[] = ":d:h:p";
        const char optstr[] = "c:t:hd:n:s:b:";
        char *ptr;

        opterr = 0;
        while ( (c = getopt(argc, argv, optstr)) != -1) {
                switch (c) {
                case 'c':
                        sprintf(check_filename, "%s", optarg);
                        printf("check_filename:%s\n", check_filename);
                        break;
                case 't':
                        sprintf(train_filename, "%s", optarg);
                        printf("train_filename:%s\n", train_filename);
                        break;
                case 'd':
			            test_dim = atoi(optarg);
                        printf("sample dim:%d\n", test_dim);
                        break;
                case 's':
			            dim_2_is_string = atoi(optarg);
                        printf("sample dim 2 is string:%d\n", dim_2_is_string);
                        break;
                case 'b':
			            boxplot_check = atoi(optarg);
                        printf("sample boxplot check:%d\n", boxplot_check);
                        break;
                case 'n':

			best_n_cluster = atoi(optarg);
                        printf("best_n_cluster:%d\n", best_n_cluster);
                        break;
                case 'h':
                        printf("Usage: ./main -t <train_filename> -c <check_filename> -d dim -s dim_2_istring -b boxplotcheck\n");
                        printf("\tsuch as:./kmeans_test -t 1000_username_new.txt -d 2 -c 1000_username_attack.txt -b 1\n");
                        printf("\tsuch as:./kmeans_test -t 1000_username_new.txt -d 2 -c 1000_username_attack.txt -b 0\n");
                        printf("\tsuch as:./kmeans_test -t 1000_test_circle.txt -d 2 -s 0\n");
                        exit(0);
                case '?':
                        printf("unknow option %c\n", optopt);
                        return -1;
                }
        }
        if (optind != argc)
                return -1;

        return 0;
}

// Makefile: gcc -g kmeans.c -DKMNEAS_TEST -lm -o kmeans_test
int main(int argc, char *argv[])
{
	int ret;
	kmeans_model_t m;	

	if (_parse_opt(argc, argv)) {
		printf("run -h for help\n");
		return -1;
	}
	if (strlen(train_filename) == 0) {
		printf("no train file, run -h for help\n");
		return -1;
	}

	memset(&m, 0x0, sizeof(kmeans_model_t));

	ret = km_read_samples(train_filename, test_dim);
	if (ret <= 0) {
		printf("no train data in %s\n", train_filename);
		return -1;
	}

	ret = kmeans_train(test_dim, km_samples, km_sample_cnt, &m, 0);
	if (ret < 0) {
		printf("train failed\n");
		goto out;
	}

    /* need to read sample again, the prev sample is preprocessed */
	ret = km_read_samples(train_filename, test_dim);
	if (ret <= 0) {
		printf("no train data in %s\n", train_filename);
		return -1;
	}
	kmeans_check(&m, km_sample_cnt, km_samples, 1, 0.0, boxplot_check);


	if (strlen(check_filename)) {
		ret = km_read_samples(check_filename, test_dim);
		if (ret <= 0) {
			printf("no test data in %s\n", train_filename);
			return -1;
		}

		kmeans_check(&m, km_sample_cnt, km_samples, 1, 0.0, boxplot_check);
	}



out:
	return ret;
}
#endif

